      SUBROUTINE INTRPL ( L, X, Y, N, U, V, IERR )
C
C  INTRPL INTERPOLATES A FUNCTION Y(X) FROM A SET OF VALUES MONOTONICALY
C  INCREASING IN X.
C
C  INPUTS :
C    L = NUMBER OF DATA POINTS (>= 2)
C    X = ARRAY OF X VALUES IN INCREASING ORDER
C    Y = ARRAY OF Y VALUES
C    N = NUMBER OF POINTS TO BE INTERPOLATED
C    U = ARRAY OF X LOCATIONS FOR OUTPUT (N VALUES)
C
C  OUTPUTS :
C    V = ARRAY OF INTERPOLATED Y VALUES
C    IERR = ERROR INDICAT0R...
C             0 = NO ERROR
C             1 = L TOO SMALL
C             2 = N LESS OR EQUAL TO 0
C             3 = X VALUES NOT MONOTONICALY INCREASING
C
      DIMENSION X(L), Y(L), U(N), V(N)
      REAL M1,M2,M3,M4,M5
      EQUIVALENCE (X2,A1,M1),(X5,A5,M5),(Q1,T3)
      EQUIVALENCE (J,SW,SA),(Y2,W2,W4,Q2),(Y5,W3,Q3)
C
      IERR = 0
      LO  = L
      LM1 = LO - 1
      LM2 = LM1 - 1
      LP1 = LO + 1
      NO  = N
      IF ( LM2 .LT. 0 ) THEN
         IERR = 1
         RETURN
      ENDIF
      IF ( NO .LE. 0 ) THEN
         IERR = 2
         RETURN
      ENDIF
      DO 11 I = 2, LO
         IF(X(I-1) .GE. X(I)) THEN
            IERR = 3
            RETURN
         ENDIF
11       CONTINUE
      IPV = 0
      DO 80 K = 1, NO
         DX = U(K)
         IF (LM2 .EQ. 0) THEN
            I = 2
         ELSE IF (DX .GE. X(LO)) THEN
            I = LP1
         ELSE IF (DX .LT. X(1)) THEN
            I = 1
         ELSE
            IMN = 2
            IMX = LO
21          I = (IMN + IMX)/2
            IF (DX .LT. X(I)) THEN
               IMX = I
            ELSE
               IMN = I + 1
            ENDIF
            IF (IMX .GT. IMN) GO TO 21
            I = IMX
         ENDIF
         IF (I .EQ. IPV) GO TO 70
         IPV = I
         J = I
         IF (J .EQ. 1)J=2
         IF (J .EQ. LP1)J=LO
         X3 = X(J-1)
         Y3 = Y(J-1)
         X4 = X(J)
         Y4 = Y(J)
         A3 = X4 - X3
         M3 = (Y4 - Y3)/A3
         IF (LM2 .EQ. 0) THEN
            M2 = M3
            M4 = M3
            GO TO 45
         ENDIF
         IF (J .NE. 2) THEN
            X2 = X(J-2)
            Y2 = Y(J-2)
            A2 = X3 - X2
            M2 = (Y3 - Y2)/A2
            IF (J .EQ. LO) THEN
               M4 = 2*M3 - M2
               GO TO 45
            ENDIF
         ENDIF
         X5 = X(J+1)
         Y5 = Y(J+1)
         A4 = X5-X4
         M4 = (Y5 - Y4)/A4
         IF (J .EQ. 2) M2 = 2*M3-M4
45       IF (J .LE. 3) THEN
            M1 = 2*M2 - M3
         ELSE
            A1 = X2 - X(J-3)
            M1 = (Y2-Y(J-3))/A1
         ENDIF
         IF (J .GE. LM1) THEN
            M5 = 2*M4 - M3
         ELSE
            A5 = X(J+2) - X5
            M5 = (Y(J+2) - Y5)/A5
         ENDIF
C
C --- NUMERICAL DIFFERENTIATION
C
         IF (I .EQ. LP1) GO TO 52
         W2 = ABS(M4-M3)
         W3 = ABS(M2-M1)
         SW = W2 + W3
         IF (SW .EQ. 0) THEN
            W2 = .5
            W3 = .5
            SW = 1.0
         ENDIF
         T3 = (W2*M2 + W3*M3)/SW
         IF (I .EQ. 1) GO TO 54
52       W3 = ABS(M5-M4)
         W4 = ABS(M3-M2)
         SW = W3 + W4
         IF (SW .EQ. 0.) THEN
            W3 = .5
            W4 = .5
            SW = 1.0
         ENDIF
         T4 = (W3*M3 + W4*M4)/SW
         IF (I .EQ. LP1) THEN
            T3 = T4
            SA = A2 + A3
            T4 = .5*(M4 + M5 + A2*(A3-A2)*(M2-M3)/SA**2)
            X3 = X4
            Y3 = Y4
            A3 = A2
            M3 = M4
         ENDIF
         GO TO 60
54       T4 = T3
         SA = A3 + A4
         T3 = .5*(M1+M2+A4*(A4-A3)*(M3-M4)/SA**2)
         X3 = X3 - A4
         Y3 = Y3 - M2*A4
         A3 = A4
         M3 = M2
C
C --- DETERMINATION OF COEFFICIENTS
C
60       Q2 = (2.0*(M3-T3)+M3-T4)/A3
         Q3 = (T3+T4-2.*M3)/A3**2
C
C --- COMPUTE POLYNOMIAL
C
70       DX = DX - X3
         V(K) = Y3+DX*(Q1+DX*(Q2+DX*Q3))
80       CONTINUE
      RETURN
      END
C
C---END INTRPL
C
